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Abstract 

We discuss the study of topological characteristics of random fields 
that are used for numerical simulation of oil and gas reservoirs and 
C/^ numerical algorithms, for computing such characteristics, for which we 

, demonstrate results of their applications. 
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^ 1 Introduction 

For the efficient extraction of oil (or gas) from oil and gas reservoirs modern 
^ technology is needed to monitor the development of fields and, in particular, 

CO the methods of geological and hydrodynamic modeling and geosteering. Now 

^ for reproduction of the real structure formation there are used probabilistic 

methods of digital geological modeling [Ij, which are as follows: 
^ An oil (gas)-bearing bed is discretized, i.e. is represented by a grid, i.e. 

^ a cover by disjoint cells which in practice often consists of parallelepipeds of 

50 m horizontal breadth and 0.4 m vertical depth. 
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Further, to each ceh, through stochastic modehng, based on concep- 
tual structure formation and statistical evaluation of the input geological 
information, there are attributed its capacitive-filtration properties such as 
porosity, permeability, compressibility, etc. 

In the development of a reservoir to control the flow of a fluid in it 
there are used hydrodynamic simulations which are software products for 
numerical solving multiphase filtration equations for which there no ana- 
lytical solutions are available. Numerical calculation is extremely resource 
consumptive and this factor is decisive in modeling with more cells, and, in 
particular, determines the size of the integrated (in the "upscaling") com- 
putational cells, as well as the need for further zeroing some of them. More- 
over, for adequate reproduction of the flow it is obviously necessary to use 
averaging filter equations and determination of effective capacitive-filtration 
characteristics of the design grid cell. 

Hence, the main purpose of a digital geological modeling in the oil indus- 
try is the creation of a geological reservoir model by regard to its geometric 
characteristics and by determining reservoir properties. Therewith it must 
be assumed that the accuracy of the resulting model is determined by the 
hydrodynamic model of the flow of a fluid in the reservoir (we want to em- 
phasize uselessness of excessive detail.) In practice the process of creation 
of a geological and hydrodynamic formation model is consecutive and often 
independent, i.e., first by geologists and then by developers. At one stage 
some design principles are used and on the other - the other, and they often 
contradict each other. Therefore, it is necessary at an early (conceptual) 
stage of a geological modeling to take into account the development data, 
and of course, they are unavoidable in the construction of a digital geological 
model. The natural question arises: how to link geology and field develop- 
ment. This is a central issue that motivates our research. Unfortunately 
there is no explicit answer to this question But it is clear that we cannot do 
without a list of characteristics of random fields (digital geological models). 

We recall some geometrical characteristics of three-dimensional digital 
geological and reservoir simulation models: the distribution of the lengths of 
streamlines [2J, the decline rate of a well discharge at a constant depression 
[3], and the fractal dimension of the model and its percolation properties 
[1]. The present work is devoted to the study of topological characteristics 
of random fields (geological and hydrodynamic models) and their relation 
to the solution of filtration equations (the data of field development). 

The problem of efficient computation of numerical characteristics de- 
scribing the topological structure of complex geometrical objects is a subject 
of computational topology that is being actively developing in recent years 
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[3 El El- In this case, often the objects under study are given as a series of 
nested into topological spaces (i.e., by filtration), and we have to trace how 
the topological structure of subspaces changes in the process of getting the 
original space. In topology, and more precisely in Morse theory, we have to 
deal with the filtration of topological spaces, which arises when considering 
the excursion sets of a function /; i.e., the sets of the form {/ > cq = const }. 
In this situation, we have to follow the change of the topological structure 
of the excursion set when cq is changed. 

In applications / by itself can have a probabilistic nature, and, so, there 
appears a problem of statistical evaluation of the impact of the excursion 
level on the topological characteristics of excursions; 

there is a problem of distinguishing topological invariants that are per- 
sistent under small perturbations of the excursion level. 

A useful tool for the study of these and other related problems are ho- 
mology and Betti numbers, and for investigation of the dependence of topo- 
logical properties of the excursion set on its level one can use persistent ho- 
mology and persistent Betti numbers. Roughly speaking, the persistent ho- 
mology estimate the portion of homology that "survives" for a given change 
of the level of the function. A detailed account of persistent homology and 
applications can be found in[8l|9l[l0l[IIl[ia[l3l[H[l5l[l6]^ 

When modeling the formation there naturally arises the permeability 
function defined by its values in each grid cell. The excursion set {/ > cq} 
of this function is a three-dimensional body modeling a reservoir for a given 
threshold of permeability. We note that this definition of a reservoir as 
the excursion set of the permeability function carries certain dangers in 
simulation, since it is difficult to clearly specify the excursion level in which 
we distinguish permeable and impervious areas, especially when we consider 
the probabilistic nature of this function. Thus, the method of comparison of 
implementations must be stable under fluctuations of excursion levels that 
lead us to use persistent topological characteristics. We note that usually 
various applications of persistent homology are due to resistance to noise at 
changing objects. 

In conclusion, we note that we study topological characteristics of real- 
izations of a random field: calculation is made after the choice of an imple- 
mentation and the excursion level. It is a reasonable problem of computing 
characteristics with taking into account the probabilistic nature of the ob- 
ject, i.e. estimation of the topological characteristics of excursions sets of 
the random field that models the structure formation. Problems of the kind 
were studied in p2|. An important problem that arises is the formulation of 
filtration equations for a random field and the relation of solutions to these 
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equations to the characteristics of the random field and that is a subject for 
further research. 



2 Computation of Betti numbers 



In [18] there is presented a numerical algorithm for computing topological 
invariants of three-dimensional bodies by using a discrete version of Morse 
theory. These invariants are the Betti numbers bo^bi, and 62, i-e. the num- 
bers of connected components, of independent one-dimensional cycles and of 
"voids" in the body. These characteristics have clear interpretations in terms 
of the permeability of a reservoir: the connectedness and compartmental- 
ization of a reservoir play primary roles in problems of its development. In 
Fig. [l]and Fig. [2] we present different realization of the same reservoir that 
are obtained by different methods cahed SGS and SPECTRAL. The SGS 
method (a successive Gauss simulation [19j) is widely accepted and is based 
on the assumption that geophysical fields are stationary both laterally and 
vertically. The method SPECTRAL was presented in [20^^ ^2T], and its main 
difference consists in the following representation of a geophysical field: 

k 

where x and y are the lateral variables, h is the vertical variable, L]^{h) are 
the Legendre polynomials, and the random processes ak(x,y) are assumed 
stationary. 

Here GL (the gamma logging) is the natural radioactivity of formation, 
and a reservoir is modeled as the excursion set {aGL < const } of the 
function 

aGL = 



GL — GL min 



GL niax GL niin 

defined on the cube of size 120 x 120 x 490. We note that we have the 
reverse inequality in the definition of an excursion because the permeability 
of a formation is inverse to its radioactivity. In Fig. [3] and Fig. |4] there 
are displayed the excursion sets {aTK < 0.6} for the two realizations given 
in Fig. [1] and Fig. [2j Different colors correspond to different connected 
components. 

For every realization the Betti numbers and the Euler characteristic 
X = 60 — ^1 + ^2 are computed. Table [l] demonstrates the dependence 
of the Betti numbers and the Euler characteristic on the excursion level. 
For every excursion level we give results of computing the topological char- 
acteristics of two different models of the same reservoir that are obtained 
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Figure 2: Realization of an oil reservoir by SGS. 



Figure 3: Excursion of a realization of oil reservoir obtained by SPECTRAL. 




Figure 4: Excursion of a realization of oil reservoir obtained by SGS. 
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by the method SPECTRAL (the upper hne) and by the SGS method (the 
lower hne). The last column contains the duration of computations on the 
processor Intel(R)Core^^i7 3.33GHz. In Fig. [sjjs] one may find graphs of 
different characteristics of excursion for both methods. We note that Betti 
numbers may distinguish the models of reservoirs obtained by the different 
methods of geostochastic modeling from the same geophysical data. 
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Table 1: The Betti numbers and the Euler characteristic. 



Remark. A computation of topological characteristics demonstrates the 
difference between the methods of geostochastical modeling, i.e. between 
SGS and SPECTRAL: the Betti numbers for different models of the same 
reservoir may vary upto 2 — 6 times. 

3 Persistent homology 

Rigorous exposition is given, for instance, in [22l [23] of cell complexes and 
the basic ideas and constructions of Morse theory that we use in the sequel. 
For computing the topological characteristics of a space it is convenient 
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Figure 7: The 2-nd Betti number. 
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to represent the space as a union of elementary "bricks," i.e. cells, which are 
"correctly" glued to each other. The resulted space is called a cell complex. 
By a 0-dimensional cell one means a point, and a union of finitely many 
0-dimensional cells forms the 0-th skeleton of a cell complex X. Let us 
consider a family of 1-dimensional cells, i.e. intervals glued to X so that the 
ends of the intervals are identified with certain 0-cells. The resulted space 
would be the 1-dimensional skeleton X^. We construct the complex X by 
successively gluing z-dimensional discs to (z — 1) -dimensional sceleta X'^~^. 

For our purposes it is enough to use cubic complexes, i.e. such cell 
complexes that all z-cells are z-dimensional cubes are glued to as fol- 

lows: every boundary face of a i-cell is an {i — l)-dimensional cube which is 
identified with some (z — 1) cube from 

Let us consider the filtration of a cell complex X by cell subcomplexes: 

= Xo C Xi C . . . C = X. 

We consider homology with coefficients in the residue group Z2. The 
filtration defines the chain of homomorphisms of the homology groups Hq = 
Hq{Xp): 

O^H^^Hl^ ...^H^ ^ H^^^ = 

for every q > 0. The compositions of successive homomorphisms from the 
chains give rise to the homomorphisms 

By definition, the persistent homology groups of dimension q are the 
groups 

Hi^j = Im/^'^' for < z < j < n + L 

Respectively by q-th persistent Betti numbers we mean the ranks of the 
persistent homology groups: bq^ = rankiJg-^. In particular, Hq^ = H^. 

Let us fix q and choose a basis . . . , e^.} for such that for 

every 1 < k < mi^ fg^^^i^k) ^ ^1^^' • • • ' ^m^+i} every I < k < nii and 
/r^'(4) = fr^\4,),k ^ k' if and only if /^''+'(e^j = 0. Hence i/^^^ 
consists of such elements that do not vanish, i.e. survive. Respectively 
the persistent homology group Hq^ consists of elements G that survive 
up to Hq. 

There is a useful graphical representation for persistent homology that is 
called a barcode ^ilOiilSj. Namely, given the dimension let us consider a 
basic element that is not an image of any element from Hq~^. Then here 



10 



exists a minimal value j > i such that fq^ (el,) = 0. Then we correspond 
to the interval (i, j). A disjoint union of all such intervals is usually 
portrayed on the two-plane by intervals parallel to the Ox axis and forms 
the g-barcode. It gives a visual representation for changing of topology of 
Xi with increasing of i. 

4 Computation of homology 

In this section we demonstrate the main ideas of the numerical algorithm for 
computing the Betti numbers of three-dimensional bodies which is presented 
in [18] by using an example of computing the persistent 0- and 2-homology 
and present some results of computations. 

Let us consider some cubic domain, in the Euclidean space, 

K= [0,7V] X [0,7V] X [0,iV] C R^ 

with some natural N. By an elementary interval / C M we mean a set of 
the form 

I = [1,1 + 1], 

where / is some natural number. Analogously we define natural square 

Q = /i X /2 c M^ 

and elementary cube 

C = /i X /2 X /3 c M^ 

where //e, /c = 1, 2, 3 are elementary intervals. Hence the domain K consists 
of elementary cubes. 

Let Ml, . . . , Mn be 3-dimensional bodies formed by elementary cubes 
and lying inside K. We assume that every Mi is the excursion set {/ > q} 
for some continuous function / defined on elementary cubes from K and 
Ci > Cj for i < j. Hence we have the filtration 

Ml C M2 C . . . C Mn. 

Variation of an excursion level from ci to results in variation of the 
topology of the excursion sets and that may be described in terms of the 
persistent homology HI = H^(Mi). 

By applying, if need be, the preprocessing of M [18j, we assume that 
two elementary cubes from M may not touch each other only at a vertex or 
along an edge. In applications that means that oil may pass from one cell 
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to another only through a common 2-dimensional face and there is no oil 
passing through common vertices and edges. 

In [18j there is proposed a numerical algorithm for computing the ho- 
mology groups of Mi by using a discrete version of Morse theory. Let us 
briefly expose the main constructions. We consider the "diagonal" linear 
function on K\ 

and the excursion sets 



Mf = {x G Mi\f{x) < a). 



A critical point of / is a vertex v ^ Mi, i.e. an integer- valued point of 
the rectangular lattice in such that when a passes af{v) the topology of 
M^ changes. All combinatorial types of critical points v = (^1,^2,^3) are 
classified in terms of their elementary neighborhoods: 

N{v) = {x G M\\xi - ki\ < l,z = 1,2,3}. 

A nondegenerate critical point has index 0, 1, or 2 being the dimension of a 
cell that glued to M^ when a passes the critical level. Moreover, there is a 
degenerate critical point, the "monkey saddle," such that two 1-dimensional 
cells are glued during passing the corresponding critical level. In Fig. [9] and 
Fig. [To] there are exposed the classical critical points: the saddle defined by 

and the "monkey saddle," defined by the 



tne equation j^x^y) — x — y 
equation f{x,y) = x^ — xy^, and also their discrete analogs. 





Figure 9: The saddle and its discrete analog. 
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Figure 10: The "monkey saddle" and its discrete analog. 

In [T8] it is constructed a chain complex 

C2(M,)^Ci(M,)^Co(M,), 

consisting of vector spaces Cq{Mi) over Z2. The basic vectors in Cq{Mi) 
correspond to the critical points of index q and, moreover, the monkey saddle 
correspond to a pair of basic vectors from Ci(M^): to each monkey saddle 
V we add a Active vertex which lies above v with respect to the level 
of / and a Active edge vv' that joins v and v'). The horizontal arrows 
denote the differentials, i.e. linear operators 82 C2{Mi) Ci{Mi) and 
d\ : Ci{Mi) Co{Mi) such that Ker^i = Im^^. We have 

S, = Im{C,+i^CJ, Z, = Ker{C,^C,_i}, H^^ZjEq 

where we assume that C3 — C_i — 0. 

The differentials are constructed explicitly [18j and for a demonstrative 
example we need only the following property: 

If is a critical point of index 1, then there is a pair of sequences of ver- 
tices L{v) — [{v^ i;^, . . . , v^)^ {v^ ^ . . . , v^)] such that they contain exactly 
two critical points of index which are v'^ and and every two consecutive 
points v^,v^_^i or v^,v^_^i are connected by a negative edge, i.e. such an 
edge that the value of / at its end is less that at its starting point. Then 
di{v) +v^. 

Our task is to construct the homomorphisms (/^^ : Cq{Mi) Cq{Mi^i) 
compatible with differentials. After that, as explained in the previous sec- 
tion, we can calculate the persistent homology and construct the barcodes 
that reflect the dynamics of change of the topological structure of a Mi as 
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i increases. Immediately we understand the arising difficulty: the natural 
inclusion Mi C M^+i induces no the natural homomorphism (f'^^ : Cq{Mi) 
Cq{MiJ^i) at the critical points, and so there are no natural homomorphisms 
of homology groups induced by the embedding Mi C M^+i. This difficulty 
is overcomed by the use of the discrete gradient flow similar to that which 
was introduced in [TS] . 

First we define the gradient descent of an arbitrary graph V formed by 
edges. Namely we assume that V is obtained by an elementary descent of 
r, if 1) (r \ r') U (r^ \ r) is the boundary (possibly without vertices) of the 
elementary face, and 2) all the vertices of \ F lie on the lower levels of / 
than all the vertices of F \ F' and F' \ F is not empty. If F' is obtained from 
F by a finite sequence of elementary descents and there is no elementary 
descent for F', then we say that F' is obtained from F by gradient descent. 

Next we assume that we have a pair of three-dimensional bodies M d N ^ 
consisting of elementary cubes. It suffices to construct homomorphisms of 
chain groups for such pairs. Let Mq and Nq be the sets of critical points, of 
index of / in M and in N . Let v G Mq. Given v ^ Nq put (fo{v) = v. 
Otherwise, there is a negative edge ei, in starting at and let vi be 
its another end. If vi ^ Nq, then we put (fo{v) = vi, and etc. We obtain 
an iterative process that results in the chain ^(v) = ('^,'^1, • • • where 
V G Mq, Vk G A^o, Vi ^ Nq foT 1 < i < k, and all edges [vi, Vi^i] are negative. 
We put (fo{v) Vk- 

Let V E Ml. Let us construct Lm{v) = [{v, v^ , . . . , v^), {v, v^ , . . . , v^)]. 
To each of the sequences from Lm{v), we add the sequence ^{V^) or $(1;+), 
respectively, and obtain a new pair of sequences [(i;, . . . , i;^, . . . , v~), (i;, . . . , 
1;+, . . . , V^)]. Let us construct the graph F consisting of the edges [v^ ^ '^z+i]? 
[v^,v^_^-^], 1 < z < /c — 1, and [i;,!;^], ['y,'^^] and let F' be obtained from F 
by the gradient descent. Obviously, the vertices v~ and v'^ G Nq and their 
constituent edges cannot down below. Therefore there exists a path 7' in F' 
which connects Vp and v'^. Let us consider all critical points wi^ . . . ^wi of 

I 

index 1 in 7' and put (fi(v) — ^wi. 

i=l 

We have (^{){di{v)) = ^o('^^+'^m) = ^p^^q ^r v ^ Mi. But di{(fi{v)) = 
I 

di{^ Wi). Since the expansions for di(wi) and for di{wi^i) have a common 

i=l 

component that is a critical point of index and the field of coefficients 
is of characteristic 2, di{(fi{v)) = v~ + v'^. Hence (po{di{v)) = di{(pi{v)), 
i.e. the differentials commute with the homomorphisms og homology groups 
induced by the embeddings. 
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The commutative diagram 



Ci(M) ^ Ci(iV) 
di i i di 

Co(M) ^ Co(iV) 

enables us to compute the persistent 0-homology. 

To compute the persistent 2-homology we have to use duahty and 
to compute the persistent 0-homology for the dual space. Namely, let us 
consider the three-dimensional body = K \ M, the complement to M, 
and the function h — —f . Clearly the critical points of indices 0, 1, and 
2 of h coincide with the critical points of indices 2, 1, and of / and so 
we may reduce the computation of 2-homology and 2-barcodes of M to the 
computation of 0-homology and 0-barcodes of M^ 

In Fig. [TT]we present the barcodes of persistent 2-homology of reservoirs 
obtained by SPECTRAL (above) and SGS (below). 
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